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Abstract 

Using first-principles molecular dynamics, we study the influence of nuclear quantum effects 
(NQEs) and nonlocal exchange-correlation density functionals (DFs) near molecular dissociation 
in liquid hydrogen. NQEs strongly influence intramolecular properties, such as bond stability, and 
are thus an essential part of the dissociation process. Moreover, by including DFs that account 
for either the self-interaction error or dispersion interactions, we find a much better description of 
molecular dissociation and metallization than previous studies based on classical protons and/or 
local or semi-local DFs. We obtain excellent agreement with experimentally measured optical prop- 
erties along pre-compressed Hugoniots, and while we still find a first-order liquid-liquid transition 
at low temperatures, transition pressures are increased by more than 100 GPa. 

PACS numbers: 64.70.Ja, 61.20.Ja, 62.50.+p, 67.90. +z 
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Hydrogen, being the most abundant element in the Universe, has a prominent role in 
planetary science. Considerable attention has thus been given to the study of its phase 
diagram at high pressure, both experimentally'^' and via first-principles (FP) simulations. 
The latter have been particularly important, in many cases being instrumental in providing 
the correct interpretation of conflicting experimental results. This is well exemplified, for 
instance, by the experimental controversy over the maximum compression along the prin- 
cipal HugoniolP'^ simulations overwhelmingly favor one with a maximum compression of 
~4.3-4.4^HIl in agreement with experiments using magnetic implosions at the Z pincri^tyj 
and converging explosive-driven shock waves^"^. Another example of the predictive capa- 
bility of FP simulations is that of a maximum in the melting line of the molecular solid 
by FP molecular dynamics (FPMD)^, subsequently confirmed by measurements^'. Un- 
fortunately though, FP methods still employ questionable approximations that could affect 
their predictions, especially when effects occur near-simultaneously, such as metallization 
and molecular dissociation. 

The emerging picture of molecular dissociation in liquid hydrogen, as suggested by FP 
simulations, is shown in Fig. [I] Below a critical temperature of ~1500-2000 K, this occurs 
through a first-order liquid-liquid phase transition (LLPT) between an insulating molecular 
liquid and a conducting atomic-like liquid. The LLPT is characterized by a discontinuous 
change in the electrical conductivity with increasing pressure combined with a small volume 
discontinuity. Above the critical temperature, the electrical conductivity and dissociation is 
then continuous with increasing pressure. While different levels of theory agree qualitatively 
regarding the LLPT^ 122 ' 25 1 , the location of it is still a matter of debate. 

Obviously, the variability in results must arise from approximations employed in the 
underlying numerical methods. The two main ones typically employed (at least in the 
case of high-pressure hydrogen) are the neglect of nuclear quantum effects (NQEs) and, in 
density-functional theory (DFT) studies, the inability to fully treat electronic correlation 
(e.g., dispersion interactions) as well as approximations to exchange; the latter typically 
resulting in a strong underestimation of band gaps^l. 

The neglect of NQEs in FP simulations is typically (but not always^^SEI]) employed, 
due to increased computational demands. Neglecting them, however, has an important 
effect on the calculated LLPT transition pressures. The use of classical protons leads to 
an overestimation, since the high-frequency vibrations of the molecular bond leads to a 
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FIG. 1: (Color online) Phase diagram of hydrogen. Diamonds (black) and downward-triangles 
(green) represent PIMD+vdW-DF2 and PIMD+PBE state-points, respectively, at which the elec- 
tronic conductivity reaches 2000 (0 cm)" 1 , which separates the insulating and conducting regimes. 
Upward-triangles (red) and squares (blue) are results from previous studies using classical protons 
in FPMD+PBE and coupled electron-ion Monte Carlo simulations^^!, respectively. The orange 
circle is the original prediction for the location of the LLPT from ScandolcF^ using FPMD within 
the local density approximation, and the turquoise diamond is the experimental measurement of 
minimum metallic conductivity by Weir, et a/P^ 

large zero-point motion (ZPM), an effect which is much smaller in the atomic phase (at 
least right at dissociation). This effect has been clearly shown by previous path-integral 
molecular dynamics (PIMD) simulations using the Perdew-Burke-Ernzerhof (PBE) density 
functional (DF)P, where the LLPT is decreased by ~40 GPa with respect to FPMD+PBE 
(i.e., classical protons). However, this would predict a transition pressure around Pt = 130 
GPa at temperature T t = 1000 K, a result not supported by current experiments^^. There 
are thus additional approximations which are likely to blame. 

The vast majority of DFs used in DFT simulations of hydrogen have been based on 
either the local density approximation or the semi-local generalized gradient approximation. 
These, however, underestimate the band gap by 1-2 eV^. This means that the metallization 
pressure, directly related to the dissociation process and thus the location of the LLPT, 
will also be underestimated. As should be clear from this discussion, neglecting NQEs 
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overestimates the transition pressure, while using a local or semi-local DF underestimates 
it. Therefore, the two errors partially cancel, providing a reasonable prediction in this 
particular instance. Such a cancellation is not always as fortunate, however, as indicated by 
the calculation of other properties near metallization, such as melting [see the Supplementary 
Material (SM)^. Therefore, rigorous simulations including NQEs with local or semi-local 
DFs must be carried out with caution. 

In this Letter, we present results from FP simulations based on PIMD to treat NQEs, but 
using nonlocal DFs in DFT. These calculations remove one of the most significant approx- 
imations made in a number of previous simulations (classical protons), while at the same 
time improve over another equally important and heretofore less-considered approximation 
(local or semi- local DFs). Such calculations allow us to study molecular dissociation in 
hydrogen with previously unattainable accuracy. 

Simulations were performed via DFT, and we focused on two nonlocal DFs. We first 
chose to use the Heyd-Scuseria-Ernzerhof (HSE) DF^, which is known to have a very small 
self-interaction error^l. We also performed simulations with the vdW-DF2 DF^HIS which 
provides a reasonable description to exchange (for a semi- local functional), but moreover 
provides an improved description of nonlocal correlation (dispersion interactions) in DFT. 
Simulations with the former were performed with VASP^ and the latter with a modified 
version of Quantum ESPRESSO (QE)P. A time-step of 8 (a. u.) 1 was used in all sim- 
ulations, and the path integrals (Pis) were discretized with a Trotter time-step no larger 
than 0.000125 K _1 . After an equilibration periods of ~0.25 ps, statistics were gathered for 
simulation times of ~ 1.5-2.0 ps, corresponding to ~ 6500-9000 time steps. A Troullier- 
Martins norm conserving pseudopotentiaP^ with a core radius of r c = 0.5 a. u. was used 
to replace the bare Coulomb-potential of hydrogen in the QE simulations; a PAW 42 was 
used in VASP. System sizes ranged from 128-432 atoms (a large number of atoms has been 
previously shown to be required for the proper description of the dissociation transition in 
liquid hydrogen at lower temperatures^). All simulations were performed at the Gamma 
point. The simulations with QE were performed with a plane- wave cutoff of 1224 eV, while 
the simulations with VASP were performed with a plane-wave cutoff of 250 eV and "Ac- 
curate" settings. Finite-temperature effects on the electrons were taken into account by 
using Fermi-Dirac smearing^. While most PIMD simulations were performed with the 
standard primitive approximation^, the simulations of 432 atoms at temperatures of 1000 
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FIG. 2: (Color online) Comparison of the measured reflectance along pre-compressed Hugoniots 



pressure-temperature curves reported by Loubeyre et al. The two sets of data correspond to 
experiments with different initial conditions. Red circles and blue crosses correspond to PIMD 
simulations with vdW-DF2 and with the optical properties being calculated using either HSE or 
PBE, respectively. The influence of the underestimation of the bandgap on optical properties 
calculated with PBE is clear. Green stars correspond to PIMD simulations with PBE and optical 
properties calculated with HSE. In this latter case, the influence on the optical properties caused 
by structural differences (from the trajectories) is far more important. 

and 1500 K utilized the accelerated PIMD method of Ceriotti, et a/.^, based on a gen- 
eralized Langevin dynamics (GLE) and the Born-Oppenheimer approximation. The use 
of the PI+GLE method was carefully tested under the relevant pressure and temperature 
conditions, in order to guarantee proper convergence^. 

As detailed further in the SA/pSl, we first compare the structure and equation of state 
(EOS) data for systems of classical protons calculated using vdW-DF2 and HSE (the latter 
more computational demanding)^. The two DFs provide pair correlation functions (PCFs) 
in very good agreement, even though pressures from HSE are ~7% smaller. This means that 
the LLPT lines predicted by the two methods will only be slightly different. This should be 
compared to the predictions of PBE, which results in a very large disagreement relative to 
either nonlocal DF. 

We performed an extended set of FPMD simulations with classical protons as well as 
PIMD simulations using vdW-DF2 at temperatures ranging from 800 to 8000 K at densities 



from Loubeyre et al.f^ to those calculated in this work. The reflectivity was calculated on the 
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in the region near molecular dissociation. At every density and temperature, optical proper- 
ties were calculated within the Kubo-Greenwood formulationpZl, by performing excited state 
calculations on 15 statistically-independent proton configurations. Note that trajectories 
and optical properties were not necessarily calculated using the same DF, the former which 
is denoted in the following notation. Results are reported in Fig. [2j in comparison to ex- 
perimental results for hydrogen along pre-compressed Hugoniots^SI. We first note that the 
reflectivity data relative to configurations obtained with PIMD-vdW-DF2 are in good agree- 
ment. The quality of the prediction is affected by the DF used in the optical calculation, 
HSE providing an excellent agreement with experiments, having a slightly lower reflectiv- 
ity than PBE. As discussed above though, this is not unexpected, due to the well known 
bandgap problem of local and semi-local DFs. On the other hand, reflectivity results from 
configurations obtained with PIMD-PBE are ~3 times larger than the experimental values, 
even when the optical calculations are performed with HSE. This effect likely derives from 
the strong tendency of PBE to favor delocalized electronic states combined with its poor 
treatment of dispersion interactions, which probably results in inaccurate proton statistical 
configurations, and thus the metallization and LLPT process altogether. 

It is important to mention that the above simulation data agrees very well with the 
SESAME EOSpSES, the latter used to convert experimental shock velocity data to pressure, 
density, and temperature. For example, our present thermodynamic data (PIMD-vdW-DF2) 
predicts a pressure only slightly higher (~3-5%) than SESAME in the relevant density range. 
Further, along the T = 5000 K isotherm, the agreement is better than 1% for pressures in 
the range of the experiments (30-60 GPa) . 

Figure [3] shows a comparison of pressure versus density along the T = 1000 K isotherm 
for both FPMD and PIMD simulations using either PBEP or vdW-DF2. Notice that both 
DFs show a plateau in the pressure, a clear indication of a first-order LLPT. There is, 
however, a further qualitative similarity in that the transition occurs between an insulating 
molecular liquid and a conductive atomic-like liquid. There is a large quantitative difference 
in the transition pressures. The inset of Fig. [3] shows a comparison of the PCF between 
FPMD and PIMD simulations using vdW-DF2. As can be seen, NQEs have a strong in- 
fluence on the properties of the molecular peak, ZPM producing a wider distribution of 
bond distances. This results in a destabilization of the molecular state, explaining the lower 
transition pressures. (Notice that the primary vdW-DF2 results shown in the figure are 
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FIG. 3: (Color online) Comparison of pressure isotherms at 1000 K for classical protons with PBE 
(red), quantum protons with PBE (green), and quantum ions with vdW-DF2 (blue). The horizontal 
dashed black lines show the pressure plateau where the LLPT takes place. Inset: Comparison of 
PCFs between simulations of classical (red) and quantum (blue) protons using vdW-DF2 at a 
density of p ~ 0.88 g/cm 3 . 

performed with PIMD, so systems of classical protons are expected to exhibit even higher 
transition pressures, above 365 GPa). 

Figure [4] shows the electronic conductivity as a function of pressure along various 
isotherms, comparing both PBE and HSE. Note that in both cases, proton configuration 
were generated with vdW-DF2. Notice also that while the conductivity values differ between 
HSE and PBE, they nonetheless agree on the existence of a jump at T = 1000 K. 

Returning to Fig. [TJ a schematic phase diagram of hydrogen in the regime of molecular 
dissociation and below T = 6000 K can be seen. The previously reported LLPT, obtained 
with classical protons and either from FPMD+PBE or coupled electron-ion Monte Carlo 
(CEIMC^ 11 calculations are shown 54 . Both vdW-DF2 (present work) and CEIMC calcu- 
lations show a considerable increase in the transition pressures with respect to PBE, with 
those from vdW-DF2 being considerably higher. Above the critical point, state-points of 
an electronic conductivity of a = 2000 (Ocm)" 1 , separating the insulating from metallic 
liquid^, are also reported using either vdW-DF2 or PBE. Loubeyre et alW^ reported that 
the metal-to-insulating threshold was located at conditions of 10% reflectivity, since accord- 
ing to the Drude model, this corresponds to an ionization of 1%. The present criterion for 
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FIG. 4: (Color online) Electronic conductivity as a function of pressure, along various isotherms. 
Results calculated with both HSE (black-lower) and PBE (red-higher) are shown for comparison. 

metallic behavior is different though. For example, from our reflectivity data, a minimum 
metallic conductivity of a = 2000 (f2cm) _1 corresponds to a reflectivity of ~0. 35-0. 40 which 
is closer to 70% of its saturation value (~0.6). This explains why our threshold line is in 
apparent disagreement with the experimental points reported in Ref. HH1 In fact, at condi- 
tions of 10% reflectivity, close to the pre-compressed Hugoniot, we observe conductivities on 
the order of a = 100-500 (ficm) -1 . Figure [l] also shows the result from the reverberation 
shock compression of S. Weir, et al^. While the temperature was not measured therein 
experimentally, but rather estimated using a model EOS, and the error bars were rather 
large, it is nonetheless clear that the presented results of the location of the LLPT and the 
dissociation regime at higher temperatures agree rather well. 

While almost all FP simulation methods agree qualitatively on the existence of a first- 
order LLPT in high-pressure hydroge n 121 1 25 1 , its precise location depends on the approxima- 
tions employed. The results reported above clearly show that NQEs and non-local DFs in 
DFT play an important role in the description of molecular dissociation and metallization. 
The two DFs considered (HSE and vdW-DF2) were originally developed with the goal of 
addressing significant limitations of local and semi-local DFs in DFT. HSE, on the one hand, 
was developed to reduce self-interaction errors in PBE in its applications to solids^. Such 
errors lead to a strong tendency to favor delocalized electronic states, which in turn lead 
to an underestimation of bandgaps by as much as 1-2 eV (in hydrogen)^. This leads to 
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a serious underestimation of the metallization pressures in both liquid and solid phases, 
and a tendency to favor metallic states (e.g., solid structures). vdW-DF and its improved 
version vdW-DF2 (employed in this work), on the other hand, were developed to account 
for nonlocal electron correlations, such as dispersion interactions in DFT. The presented 
results indicate that, at least close to dissociation, both HSE and vdW-DF2 DFs produce 
very similar structures in liquid hydrogen. Since the physical effects addressed by both DFs 
are not directly related to each other, and that both effects are expected to be relevant in the 
molecular phase, it is important to recognize that the LLPT pressures might still change if 
a DF which combines both hybrid exchange and non-local correlation were to be employed. 

The goal of this Letter was not to predict which functional (HSE and vdW-DF2, etc.) is 
more accurate, since answering that question requires the use of more accurate methods^. 
We can however mention several possibilities that explain the observed behavior, the reason- 
able agreement between either nonlocal DF as well as their large disagreements with PBE. 
Both DFs predict shorter molecular bonds compared to PBE; in the limit of low density, 
the bond length predicted by vdW-DF2 agrees very well with measured values while that 
of PBE is overestimated by ~39cP^. This is obviously an important factor on dissociation. 
Second, the exchange portion of the vdW-DF2 functional was constructed to reproduce 
exact-exchange results^, which may explain its similarity to HSE. Finally, both dispersion 
interactions and a reduced self-interaction will lead to a more stable molecular state. An 
even more promising alternative to DFT is the use of quantum Monte Carlo first-principles 
methods, for example CEIMCP} using accurate trial wave functions, such as those con- 
structed from HSE orbitals. We must also recognize that, while the use of the vdW-DF2 
DF made large improvements in the description of molecular dissociation in hydrogen near 
the LLPT, standard semi-local DFs like PBE have been shown to be successful in describing 
other materials when combined with the HSE DF for the calculation of optical properties^. 

While the presented results are obviously important to high-pressure hydrogen, they also 
suggest that NQEs will have a strong influence on the bonding properties of other hydrogen- 
rich materials, particularly in the description of transitions between phases with different 
bonding characteristics. Although this has been demonstrated in some cases, such as high- 
pressure ice^, and similar effects could be present in many other materials (e.g., methane 
and ammonia), the simultaneous and proper inclusion of NQEs has been largely ignored 
in the field of FP simulations. With the development of efficient PI methods, such as the 
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PI+GLEpH, and faster computers, we expect that future simulations will routinely include 
them, neither them neglecting or employing approximations at the harmonic level. 

In conclusion, we have studied liquid hydrogen at high pressure using nonlocal exchange- 
correlation DFs in DFT, namely HSE and vdW-DF2. Both produce similar descriptions of 
the liquid, with large increases in the LLPT pressures (where molecular dissociation occurs) 
by more than 100 GPa (at lower temperatures), from earlier studies using PBE. We also 
presented a detailed study of the influence of NQEs in the LLPT in combination with vdW- 
DF2. Remarkable agreement with experiment was observed for optical properties along 
pre-compressed Hugoniots, as well as with reverberating shock compressed measurements at 
low temperatures. The improved description further confirms the existence of a first-order 
LLPT between an insulating molecular liquid and a conductive atomic-like state at high 
pressures and below a critical temperature of T c w 1500-1000 K. Since this work presents 
a highly-accurate prediction of the location of the LLPT in hydrogen, it should serve as a 
clear goal for future experimental and theoretical works in this field. 
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Notice that HSE simulations are very sensitive to the choice of k-points and the size of the 

simulation cell, which means that the transition pressures will depend on the particular details 

of the simulation. This is not the case for the other DFs employed in this work. While the 

precise location of the transition using HSE is a matter of future work, away from the transition 

the dependence on simulation details should be considerably smaller. 

The CEIMC method is an alternative FP simulation method based on a quantum Monte Carlo 
treatment of the electronic potential energy surface, combined with a classical or quantum 
treatment of the ions using traditional Monte Carlo methods. 

Preliminary results in solid molecular hydrogen show that PIMD+vdW-DF2 is able to reproduce 
the measured bandgap of the solid, while PIMD+HSE predicts gaps that are ~1 eV too low-^M 
These results provide evidence in favor of the higher transition pressures predicted by the 
PIMD+vdW-DF2 simulations. 
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